function [x,y,u]=roms_zslice( rslice_obj, varargin)
% ROMS_ZSLICE:  Returns horizontal slice from ROMS Netcdf file.
%
% Returns a matrix containing a horizontal slice at a specified depth
% at a given time step from a SCRUM NetCDF file.  Regions of grid that
% are shallower than requested value are returned as NaNs.
%
% USAGE: [x,y,u]=roms_zslice ( rslice_obj, timestep, depth)
% USAGE: [x,y,u]=roms_zslice ( rslice_obj, varname, timestep, depth)
%    rslice_obj:     
%        global structure that contains a lot of needed information
%    varname:
%        If provided, then don't use the variable named in rslice_obj
%    time:    
%        time index, must be zero-based
%    depth:   
%        depth in meters.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% $Name: rslice-0_9_2 $
% $Id: roms_zslice.m,v 1.6 2005/08/25 15:56:57 jevans Exp $
% AUTHOR:  johnevans@acm.org
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin == 3
	varname = rslice_obj.varname;
	timestep = varargin{1};
	depth = varargin{2};
else
	varname = varargin{1};
	timestep = varargin{2};
	depth = varargin{3};
end

ncfile =  rslice_obj.ncfile_list{rslice_obj.ncfile_index};
varinfo = nc_getvarinfo ( ncfile, varname );


%
% Default return values.
x = []; y = []; u = [];

%if (nargin~=3)
%  help roms_zslice; return
%end

depth = -abs(depth); 			% account for depth or z







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get variables for z eqn [z = zeta * (1 + s) + hc*s + (h - hc)*C(s)]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% w is defined at different locations than the others
switch ( varinfo.Dimension{2} )
case 's_rho'
	sc = rslice_obj.rslice_grid.roms_grid.sc_r;
case 's_w'
	sc = rslice_obj.rslice_grid.roms_grid.sc_w;
otherwise
	fprintf ( 2, '%s:  unhandled z dimension %s.\n', mfilename, varinfo.Dimension{2} );
end


x = rslice_obj.rslice_grid.x.data;
y = rslice_obj.rslice_grid.y.data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Retrieve Zeta.  
zeta = nc_varget ( ncfile, 'zeta', [timestep 0 0],[1 -1 -1] );

hc = rslice_obj.rslice_grid.roms_grid.hc;
h = rslice_obj.rslice_grid.roms_grid.h;
Cs_r = rslice_obj.rslice_grid.roms_grid.Cs_r;

%
% Guessing at this.
if length(sc) > length(Cs_r)
	n = length(sc);
	Cs_w = zeros(n,1);
	Cs_w(1) = -1;
	for j = 2:n-1
		Cs_w(j) = ( Cs_r(j) - Cs_r(j-1) )/2 + Cs_r(j-1);
	end

	Cs = Cs_w;
else
	Cs = Cs_r;
end


[eta_rho_length, xi_rho_length] = size(x);
s_rho_length = varinfo.Size(2);
y_length = varinfo.Size(3);
x_length = varinfo.Size(4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the depth.
n = length(sc);
z = zeros(n, eta_rho_length, xi_rho_length );
for i = 1:n
	z(i,:,:) =  zeta * (1+sc(i)) + hc*sc(i) + (h - hc)*Cs(i);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reshape the depth such that each vertical profile is a column.
y_dim_name = varinfo.Dimension{3};
x_dim_name = varinfo.Dimension{4};
if ( strcmp(y_dim_name,'eta_rho') & strcmp(x_dim_name,'xi_rho') )
	;
elseif ( strcmp(y_dim_name,'eta_u') & strcmp(x_dim_name,'xi_u') )
	y = (y(:,1:xi_rho_length-1) + y(:,2:xi_rho_length))/2;
	x = (x(:,1:xi_rho_length-1) + x(:,2:xi_rho_length))/2;
	z = (z(:,:,1:xi_rho_length-1) + z(:,:,2:xi_rho_length))/2;
	zeta = (zeta(:,1:xi_rho_length-1) + zeta(:,2:xi_rho_length))/2;
	h = (h(:,1:xi_rho_length-1) + h(:,2:xi_rho_length))/2;
elseif ( strcmp(y_dim_name,'eta_v') & strcmp(x_dim_name,'xi_v') )
	y = (y(1:eta_rho_length-1,:) + y(2:eta_rho_length,:))/2;
	x = (x(1:eta_rho_length-1,:) + x(2:eta_rho_length,:))/2;
	z = (z(:,1:eta_rho_length-1,:) + z(:,2:eta_rho_length,:))/2;
	zeta = (zeta(1:eta_rho_length-1,:) + zeta(2:eta_rho_length,:))/2;
	h = (h(1:eta_rho_length-1,:) + h(2:eta_rho_length,:))/2;
else
	disp(sprintf('what are the dimensions??, %s and %s??', y_dim_name, x_dim_name ) );
	help roms_zslice;
	return;
end
z = reshape ( z, [s_rho_length y_length*x_length] );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Retrieve the requested variable at the given timestep.
data = nc_varget ( ncfile, varname, [timestep 0 0 0],[1 -1 -1 -1] );


%data = permute ( data, [3 2 1] );
data = reshape(data,[s_rho_length y_length*x_length]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define z and data at zeta (elevation)
z = [z; reshape(zeta, [1 y_length*x_length] ) ];
data = [data; data(s_rho_length,:)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define z at values greater than zeta to be Inf
% define data to be NaN
z = [z; Inf*ones(1,y_length*x_length)];
data = [data; NaN * ones(1,y_length*x_length)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define data at bathymetry values to be same as data at lowest 
% existing values.
z = [ -1*reshape(h,[1 y_length*x_length]); z];
data = [ data(1,:); data ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define z at depth greater than bathymetry to be -Inf
% define data at depth greater than bathymetry to be NaN
z = [ -Inf*ones(1,y_length*x_length); z ];
data = [ NaN * ones(1,y_length*x_length); data ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the indices of data values that have just greater depth than
% depth.
zgreater = ( z < depth );
zg_ind = diff(zgreater);
zg_ind = find(zg_ind~=0);
zg_ind = zg_ind + [0:1:length(zg_ind)-1]';
data_greater_z = data(zg_ind);
depth_greater_z = z(zg_ind);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the indices of the data values that have just lesser depth
% than depth.
zlesser = ( z > depth );
zl_ind = diff(zlesser);
zl_ind = find(zl_ind~=0);
zl_ind = zl_ind + [1:1:length(zg_ind)]';
data_lesser_z = data(zl_ind);
depth_lesser_z = z(zl_ind);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interpolate between the data values.
alpha = (depth - depth_greater_z) ./ ( depth_lesser_z - depth_greater_z );

data_at_depth = (data_lesser_z .* alpha) + (data_greater_z .* (1-alpha));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now reshape the data into a square grid.
u = reshape ( data_at_depth, [y_length x_length] );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If there is a suitable mask variable, we will want to mask out certain cells.
% For now, mask out the land cells.
if ( strcmp(y_dim_name,'eta_rho') & strcmp(x_dim_name,'xi_rho') )
  %Get the mask_rho variable

    mask_inds = find(rslice_obj.rslice_grid.roms_grid.mask_rho==0);
    u(mask_inds) = NaN * ones(size(mask_inds));
elseif ( strcmp(y_dim_name,'eta_u') & strcmp(x_dim_name,'xi_u') )
  %
  %Get the mask_u variable
    mask_inds = find(rslice_obj.rslice_grid.roms_grid.mask_u==0);
    u(mask_inds) = NaN * ones(size(mask_inds));
elseif ( strcmp(y_dim_name,'eta_v') & strcmp(x_dim_name,'xi_v') )
  %Get the mask_u variable
    mask_inds = find(rslice_obj.rslice_grid.roms_grid.mask_v==0);
    u(mask_inds) = NaN * ones(size(mask_inds));
else
  disp(sprintf('what are the dimensions??, %s and %s??', y_dim_name, x_dim_name ) );
  help roms_zslice;
  return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% all done with netcdf file operations
return
